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Hyperluminous infrared galaxies (HyLIRGs) are the rarest and most 
extreme starbursts and found only in the distant Universe (z 2 1). They 
have intrinsic infrared (IR) luminosities Lir 2 10” L and are commonly 


found to be major mergers. Recently, the Planck All-Sky Survey to Analyze 
Gravitationally-lensed Extreme Starbursts project (PASSAGES) searched 
~10* deg’ of the sky and found ~20 HyLIRGs. We describe a detailed study 

of PJO116-24, the brightest (UL, = 2.6 x10" Lo, magnified with p ~ 17) 
Einstein-ring HyLIRG in the southern sky, at z= 2.125, with observations 
from the near-IR integral-field spectrograph VLT/ERIS and the submillimetre 
interferometer ALMA. We detected Ha, HB, [N 11] and [S 11] lines and 
obtained an extreme Balmer decrement (Ha/Hf = 8.73 + 1.14). We modelled 
the molecular-gas and ionized-gas kinematics with CO(3-2) and Ha data 

at -100-300 pc and (sub)kiloparsec delensed scales, respectively, finding 
consistent regular rotation. We found PJ0116-24 to be highly rotationally 
supported (V,ot/0o0, mot. gas ~ 9-4) with a richer gaseous substructure than 

other known HyLIRGs. Our results imply that PJO116-24 is an intrinsically 
massive (Mbaryon = 10™° Mo) and rare starbursty disk (star-formation rate, 
SFR=1,490 M, yr”) probably undergoing secular evolution. This indicates 
that the maximal SFR (21,000 M, yr’) predicted by simulations could occur 
during a galaxy’s secular evolution, away from major mergers. 


All-sky surveys conducted with infrared (IR) satellites, such as the Infrared 
Astronomical Satellite (IRAS), the Wide-field Infrared Survey Explorer 
(WISE) and the Planck satellite, have led to great breakthroughs in finding 
rare hyperluminous infrared galaxies (HyLIRGs)’ °. The brightest HyLIRGs 
have IR luminosities 4L > 10" Land are often found to be gravitationally 
lensed. High-resolution imaging with the Hubble Space Telescope (HST) 
inthe optical and near-IR and with (sub) millimetre and radio interferom- 
eters like the Very Large Array (VLA), the Northern Extended Millimeter 
Array (NOEMA) and the Atacama Large Millimeter/Submillimeter Array 
(ALMA) have revealed that HyLIRGs are almost exclusively gas-rich merg- 
ers’ ?. In most cases, they also host an active galactic nucleus (AGN). 

A widely accepted scenario is that HyLIRGs are the distant 
higher-luminosity tail of the local ultra-luminous IR galaxies (ULIRGs; 


Lr> 10” Lo) with extreme starburst activities triggered by major merg- 
ers. Thus, they fit into the well-known merger-driven ULIRG and quasar 
paradigm”. Theoretically, however, a second possibility has been 
raised that HyLIRGs are very young, ‘primaeval’ galaxies that are going 
through their maximal star formation periods. Hydrodynamical 
simulations’*”’ of massive, turbulent, gas-rich galaxies have provided 
evidence that extreme star-formation rates (SFR = 1,000 M, yr”) can 
be achieved in the brightest submillimetre galaxies at redshift 2-3. Star 
formation is fuelled by cold gas rapidly inspiralling towards the centre 
during their secular evolution, that is, away from any major merger. 
The theories are still not robustly constrained by observations due to 
the rarity of HyLIRGs with high-resolution, high-quality, multi-tracer 
observations. 
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Table 1| Physical properties of PJO16-24 


Properties Values Units 
(Apparent, magnified by lensing) 
0.35 
log(uM.) 12.36037 (Mo) 
log(uLe) 14.4170.55 (La) 
USFRęg 25,120+6,500 (Ma yr”) (Chabrier IMF) 
log(uMisu) 12.494002 a (Mo) 
EBEV 0.50+0.13 (mag) 
E(B-V)neb, seo 1.00+0.26 (mag) 
(Intrinsic, corrected for lensing’) 
log(M.) 1116+0.35 (Mə) 
log(Lir) 13.17+0.29 (Lo) 
SFRig 1,480+400 (Mg yt") (Chabrier IMF) 
SFRw 8+2 (Moy) (Chabrier IMF) 
SFRuv» ir 1,490+400 (Mo yr”) (Chabrier IMF) 
log(Mism) 11.26+0.32° (Mz) 
(Kinematic modelling) 
log(Moisk+bulge, dyn) 11.2970? (Mo) 
+0.07 
BiT 0.09 0.01 (Mo) 
log(Mpmdyn) BBA (Mo) 
Rett 3.6700 (kpc) 
Vrot, mot. gas(Rett) 40175) (kms”) 
Oo mol. gas 42.851 (kms”) 
Vrot/ Oo, mot. gas(Rett) 9.4+0.9 
Vrot, ion. gas(Rett) 389+31 (kms") 
Oo, ion. gas 68.6477 (kms”) 
Vrot/ Oo, ion. gas(Rett) 5.7 + 0.6 
inati i 8.5 
Inclination angle i Site (deg) 
iti +0.8 
Position angle (PA) 80.4708 (deg) 


Derived from an updated run of the radiative transfer model in Harrington et al.” with 
additional CO and [C |] lines and dust continua from new observations (‘Interstellar medium 
mass from CO and [C |] line modelling’; K. Harrington et al., manuscript in preparation). "We 
corrected the stellar mass with u=16 and IR luminosity and ISM mass with u =17 based on our 
lens model and the distributions of stellar, CO and dust emission. 


A unique HyLIRG under the microscope of strong 
lensing 

In recent years, hundreds of apparent HyLIRGs have been discovered 
from =10* deg? of Planck, Herschel and South Pole Telescope (SPT) 
surveys, including PASSAGES, which searched -11,500 deg? of sky*°. 
Among these, there are only a few tens of HyLIRGs with pl, > 10* Lo, 
and those that are lensed into Einstein rings (thus, with large magnifica- 
tions) are even rarer. Here we select the brightest Einstein-ring HyLIRG 
in the southern sky from PASSAGES, PJ011646.8-243702 (hereafter 
PJ0116-24), at z= 2.125 and with 4L = 2.5 x 10 L, (see physical proper- 
ties in Table 1), for a pioneering and comprehensive study of its heavy 
dust attenuation and kinematics, about which there is presently little 
information for the HyLIRG population. We obtained high-resolution 
CO(3-2) observations with ALMA to trace the cold-gas distribution 
and kinematics at ~0.16” angular resolution and high-quality near-IR 
integral field unit (IFU) observations with a new versatile instrument, 
the enhanced resolution imager and spectrograph (ERIS)”, installed on 
the Very Large Telescope (VLT)’s UT4, to trace the ionized-gas distribu- 
tion and kinematics at ~0.68” angular resolution. We also measured the 


Table 2 | Line fluxes, Balmer decrement, dust attenuation 
and SFRs derived from ERIS IFU data 


Line name and Value Units and notes 


derived properties 


(Apparent, magnified by lensing) 


Ha 517x10 "+1.68x10" (ergs 'cm”) 

HB 5.92x10%+7.51x10" (ergs 'cm”) 

Hanae 513x10 +1.66x10"® (ergs cm”) 

HB etuncor” 5.59x108+710x10" (erg s” cm?) 

[O m] A4959 7.43x10"+4.56x10" = (ergs 'cm”) (tied to [O 111] 

5007) 

[O m] ASOO7 3.31x10+5,78x10 = (ergs cm”) 

[N 1] A6548 1.03x10"+1.53x10"® (ergs cm”) 
N 1] A6584 3.32x105+1.74x10"®  (ergs”cm?) 

[S 1] A6717 8.10x10°°+5.45x10” (ergs cm”) (tied to [S 1] A6731) 
S 11] 46731 5.73x10"8+3.85x10" ~~ (ergs cm”) 

Ha/HB 8.73 +114 

Ha/HB sunon? 919+1.20 

EW esses 0.95+013 (mag) 

Ayib 2.95+0.13° (mag) 

LUSFR io, notcom 980+30° (M, yr“) (Chabrier IMF) 

USFRija, con® 7,900+1,000° (Moyr’) (Chabrier IMF) 


(Intrinsic, corrected for lensing) 


SFRita, not corr 58+2° 


SFRita, corr 


(M, yr") (Chabrier IMF) 
(Mg yr’) (Chabrier IMF) 


Directly measured values, uncorrected for the stellar Balmer absorption. "Converted from 
Ha to SFR using the Kennicutt and Evans conversion® and an assumed Calzetti et al.” 
attenuation law with Ry=3.1 (ref. 32). “SFR inferred from the attenuation-corrected Ha line 
luminosity using the Balmer decrement. 


470+60° 


global Balmer decrement and metallicity using the observed rest-frame 
optical spectral lines (Table 2). 

Asshownin Fig. la, PJO116-24 is well resolved into a nearly complete 
Einstein ring in the HST 1.6 ppm image (ref. 6 and J. Lowenthal et al., 
manuscript in preparation), which provides rich structural informa- 
tion: two bright emission peaks A1 and A2; several stretched arcs; and 
four clump images B1 to B4. The distribution of cold gas, as traced by 
the CO(3-2) line-integrated emission (Methods), exhibits extended 
structures that surround the HST emission peaks. 

The CO(3-2) line-centre velocity map, as shown in Fig. 1b, clearly 
exhibits an ordered rotation with line-of-sight velocities ranging from 
about -330 to +330 km s” relative to the systematic redshift veloc- 
ity of PJO116-24. We performed detailed kinematic modelling and 
quantified the rotation-to-dispersion ratio aS Vrot/0o, mol. gas = 9-4 + 0.9 
for the molecular gas and V,ot/0o0, ion. gas = 5-7 + 0.6 for the ionized gas 
at the baryonic-mass effective radius of R. = 3.6 kpc (corrected for 
beam-smearing, line spread function and inclination; Methods). 
This rapid and ordered rotation is rarely found in HyLIRGs and is best 
constrained among all known HyLIRGs. The HST emission peaks, A1 
and A2, are at the near-zero velocity positions and map to the same 
source-plane (intrinsic) position, confirming that they are both the 
photometric and kinematic centres of the galaxy. 

The ionized-gas distribution, as traced by the Ha line-integrated 
emission from the ERIS IFU data (Methods) in Fig. Ic, further shows 
an interesting configuration. The gas extends across the disk, tends 
to be more pronounced on the blueshifted side of the galaxy and has 
consistent detections at all BI-B4 locations. Given that Ha emis- 
sion is sensitive to dust attenuation, the blue/red-side asymmetry 
could be the combined effect of spatially varying dust attenuation, 
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a PJO116-24—an all-sky selected, lensed, intrinsic 


hyperluminous-infrared (HyLIRG) Einstein ring 
B4 


Fig. 1| Stellar, cold-gas and ionized-gas emission distributions of the rare 
HyLIRG Einstein ring, PJ0116-24. a, Lensed distribution of the stellar light as 
traced by HST 1.6 pm emission to a Sauron colour scale with blue-white-green- 
yellow. The foreground lens has been subtracted, and the PSF = 0.2” is indicated 
at bottom left. Also shown in red is the lensed distribution of the cold gas as 
traced by the ALMA CO(3-2) line intensity. The beam, 0.19” x 0.16”, is indicated at 
bottom right. b, CO(3-2) line-centre velocity map (-330 to +330 kms from dark 
blue to light pink). c, Comparison of the distributions of the cold gas (red) and the 
ionized gas as traced by the ERIS/VLT Ha line intensity (white). d, Source-plane 
distributions (delensed from the full Einstein ring) of emission from stars, cold 


lonized gas 


© 
R 
3 


kpc = 
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gas and ionized gas. The AGN at the galaxy centre is lensed into the positions A1 
and A2as labelled ina and b. The clump at the north-east -4 kpc from the centre is 
lensed to the positions B1 to B4 as labelled in a and c. (The further north Ha blob is 
an artefact of delensing due to the noise peaks, for example, around B2, scattered 
across the caustic curves.) Caustic curves from our lens model are indicated by 
the yellow lines in d, and white lines mark the CO arm-like substructures. The 
CO(3-2) ‘delensed beams’ at the galaxy centre reconstructed from positions A1 
and A2 are indicated by the red ellipses at bottom left (FWHM = 0.09” x 0.06” 

and -0.15” x 0.04”, respectively; the delensed beam is even smaller closer to the 
caustic lines). 


the underlying distribution of star formation across the galaxy 
and lensing. 

Through detailed new lens modelling using high-resolution HST 
and CO channel maps and mesh-grid-based delensing (Methods), we 
reconstructed the source-plane CO(3-2) and Ha line-integrated emis- 
sion distributions as shown in Fig. 1d. The magnification of the total 
flux of the resolved HST and CO(3-2) emission is estimated as ~16 and 
~17, respectively, using the observed and delensed maps. Due to the 
lensing effect, the spatial resolution (beam or point spread function 
(PSF)), which is uniform in the observed maps, varies with location 
in the source-plane delensed maps. For instance, when delensing a 
Gaussian beam at the A2image of the galaxy centre, the delensed beam 
has an approximately Gaussian shape with a minor-axis full width at 
half-maximum (FWHM) of -300 pc and an axis ratio of -0.3, elongated 
along north-east to south-west. This delensed beam varies across the 
source plane. A maximum resolution of <100 pc was achieved along the 
Einstein ring when crossing the critical curve (for example, between 
A1-B4 and A2-B2; Methods). 

Under the microscope of strong lensing, we found that PJO116-24 
has the intrinsic characteristics of massive star-forming disk galaxies 
at z= 2 that have been well studied with adaptive optics at kiloparsec 
scales”. That is, it has a central stellar concentration, an extended, 
massive, cold-gas-rich disk, and giant (kiloparsec scale) star-forming 
clump(s) that formed due to the strong Toomre instability of a gas-rich 
disk. The CO-traced gas structures are within about one effective 
radius of -3.6 kpc and havea Toomre’s Q parameter Q,,, ~ 0.3 (based on 
equation (3) of ref. 19, Qas = A / faas X (Oo / Vrot), Where a = 1.4 is a factor for 
ageneral flat rotation curve, f,,, = 0.56 is the gas fraction and V,o./0 mol. gas 
= 9.4 is from our kinematic modelling). This Q,,, is well below the typi- 
cal stability criterion for a thick disk Qor = 0.67 (ref. 19). Thus, the gas 
naturally fragments into structures at (sub) Toomre scales or (sub) 
kiloparsec sizes (A; ~ fgasRaisx: ref. 20). In our work, we resolved very 
rich substructures down to ~100-300 pc within the most rotationally 
supported HyLIRG. 


We performed rest-frame ultraviolet (UV) to millimetre spectral 
energy modelling (Methods) and found that PJ0116-24 is forming stars 
at arate of SFRuv+r = 1,490 + 400 M, yr”, rivalling the small number of 
known highest-SFR systems (-1,000-3,000 M, yr”) at high redshift” ™*. 
The non-major-merger scenario of PJ0116-24 is further supported by 
the cold-gas distribution, which exhibits a deficit within the central 
~500 pc but is enriched at 1-2 kpc and extends out to ~3-4 kpc with 
stream- or arm-like features (as outlined by the white lines in Fig. 1d). 
These features can naturally form in turbulent, gas-rich, secularly 
evolving disks when gas is accreted and transported into the centre”. 

Agiant clump at ~4 kpcin projection north-east of the galaxy centre 
inthe source plane is lensed into the BI-B4 emission peaks. It is detected 
in all stellar, Ha and CO tracers and is probably a kiloparsec-scale, 
young, star-forming clump with either stellar or gas mass accounting 
for only <3% of that of PJO116-24 (Methods). With FWHM = 1.2 kpc for 
a two-dimensional (2D) Gaussian fit of the delensed HST image, it is 
x10 to x100 larger than the star-forming complexes in the Milky Way 
but typical of giant clumps in massive star-forming galaxies at z = 2-3. 
These clumps also naturally form due to the galaxies’ high gas frac- 
tions (fzas = Myas/(Myas + M.) 2 0.5)”° and increased Toomre instability 
(refs. 19,20,26-29; see also simulations in refs. 25,30,31 and Methods). 

Figure 2 provides a clearer view of the delensed distributions of 
cold andionized gas. The two gas tracers exhibit remarkably consistent 
velocity, velocity dispersion and position-velocity (P-V) maps, albeit 
apparently affected by corresponding beam-smearing, which thus 
needs detailed dynamical modelling. We performed state-of-the-art 
dynamical modelling based on direct fitting of the high-resolution 
CO data in the image plane (Methods), which reproduced well the P-V 
diagrams of both CO and Ha under different beam-smearing condi- 
tions (Fig. 2g,h). The centrally peaked apparent velocity dispersion is 
also strong evidence for an ordered disk rather than a major merger". 
Additionally, the clump appears to follow the systematic rotation and 
has an indistinguishable cold-gas velocity dispersion as the disk, again 
inline with anon-merger scenario. 
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Fig. 2 | CO and Ha line intensity and kinematic maps of PJ0116-24 after 
correcting for lensing. a,b Line intensity maps for CO (a) and Ha (b). 

c,d, Velocity maps for CO (c) and Ha (d). e,f, Velocity dispersion maps (vel. disp. 
maps) for CO (e) and Ha (f). g,h, P-V diagrams for CO (g) and Ha (h). HST stellar 
light distribution is indicated by the contours in a-f. The angular resolution for 
reconstructions from images Aland A2 is indicated by the ellipses at the lower 
left corners of a and b. Colour bar next toc and d displays the velocities in units 
of km s”. Colour bar next to e and f displays the velocity dispersion in units of 
km s™. In gand h, contours represent our best-fitting kinematic model’s P-V 
diagrams with the same beam-smearing effect as the data. Dotted lines indicate 
our best-fitting intrinsic rotation curves without beam-smearing but with the 
inclination effect. Thus, they match the contours where the rotation curve 
flattens and the beam-smearing is minimized. Dashed lines indicate the intrinsic 
rotation curve after beam-smearing and inclination corrections. In h, the [N 11] 
6584, Ha and [N 11] 6548 lines (from top to bottom) are shown in the same map 
due to their proximity in wavelength and velocity. The vertical axis of hindicates 
the velocity with zero centred at the Ha line. 


Balmer decrement, metallicity and ionization 
source 

Our ERIS Science Verification Program (principal investigator (PI) 
D. Liu) observed PJO116-24 in two near-IR bands: the H-band cover- 
ing the redshifted HB and [O 111] AA4959,5007 lines and the K-band 
covering the Ha, [N 11] AA6548,6584 and [S 11] AA6717,6731 lines. 
All lines were detected with signal-to-noise ratio (S/N) = 2-20. 
This dataset provides a full suite of strong-line diagnostics for 
studying the galaxy’s dust attenuation, metallicity and ionization 
source. 


The ratio of Hato HB line fluxes, known as the Balmer decrement, 
isa reliable indicator of the dust attenuation of the emitting source. 
In case b recombination with an electron temperature 7, = 10+ K 
and density n, = 10?-10* cm”, appropriate for star-forming regions, 
the intrinsic, unattenuated Ha/Hf flux ratio is 2.86. In compari- 
son, we observed an attenuated Ha/Hf flux ratio of 8.73 + 1.14 (cor- 
rected for stellar Balmer absorption). Assuming a Calzetti et al.” 
attenuation curve with R = 3.1 for the nebular line, we obtained 
a nebular reddening of A), neb = 2.95 + 0.13 mag. This is consistent 
with the dust attenuation determined through our fitting of the 
spectral energy distribution (SED; Methods). Converting the Ha 
line flux to an intrinsic, attenuation-corrected SFR”, we obtained 
SFRua, corr = 470 + 60 M, yr”, -31% + 9% of the total SFRuv+r. The appar- 
ent, uncorrected SFRya, notcorr ~ 58 + 2 M, yr ‘is as low as ~3.9% + 1.1%. 
The small fractions indicate that a substantial amount of Ha emission 
is too obscured by dust to be detected. The high magnification in this 
Einstein ring helps to unveil the less obscured parts of the galaxy, but 
the measured Balmer decrement places only a lower limit on the total 
dust obscuration for the whole galaxy. 

Figure 3 shows the ERIS H-band HB, [O 111] doublet, and K-band 
Ha and [N 11] doublet spectra, as well as a Balmer decrement versus 
stellar mass diagram. PJO116-24 occupies a unique parameter space 
distinct from the vast majority of local galaxies in the Sloan Digital Sky 
Survey (SDSS) and z= 2 main-sequence star-forming galaxies in the 
MOSDEF survey (see Methods for detailed comparisons). The lack of 
galaxies in this unique parameter space of high M. and high Ha/Hf is 
for several reasons: the decline of stellar mass function at the massive 
end, the decline of the IR luminosity function at the bright end and the 
challenge of performing rest-frame optical spectroscopy for heavily 
obscured galaxies like PJO116-24. The PASSAGES sample will be ideal 
to fill this parameter space thanks to the strong lensing that makes 
rest-frame optical IFU imaging spectroscopy possible. Figure 4 further 
presents several sets of metallicity diagnostics using various line ratios 
as labelled in each panel. Using the Curti et al.** calibrations, which 
are shown as the solid blue lines in Fig. 4, we interpolated the metal- 
licity 12 + log(O/H) from the observed line ratios. These diagnostics 
indicate solar to supersolar metallicity. This is much higher than in 
non-starburst galaxies at the same redshifts, which is typically ~1/10 
solar”. The N2 (Ha/[N1I]) and RS32 ([S 11]/HB + [O 111]/Ha) diagnostics 
indicate an ~0.2 dex supersolar metallicity. Whether the supersolar 
metallicity from the only two calibrations is reliable or not is still an 
open question awaiting more HyLIRGs studied in the same way. 

Our analysis using the BPT diagnostic diagram (as detailed in 
Methods) further reveals that PJO116-24 is dominated by star forma- 
tion, with alow contribution from the AGN, in line with recent studies 
using X-ray spectroscopy” and high-resolution radio observations‘. 
For the highest-S/N Ha and [N 11] lines, their line profiles are spec- 
troscopically well resolved, as shown in Fig. 3, and our analysis indi- 
cates that the complex line profile composed of anarrow component 
(FWHM = 240 km s”) and a broad component (FWHM = 510 km s”) is 
mainly due to the asymmetric lensing of the rotating disk (with the blue 
side being more magnified; Methods). 


Central exhaustion of cold gas and the growth of 
the bulge 

The centre of PJO116-24 exhibits characteristics of negative ‘feedback’. 
That is, star formation is suppressed due to the lack of cold gas within 
the central half kiloparsec, which could be attributed to past activities 
of an AGN?” or acircumnuclear starburst”, The central stellar 
mass surface density of PJO116-24 is X., <0.5kpc % 3 X 10°° Mo kpc 
within -0.5 kpc, or ©, <ıkpe & 15 X 10°° Mo kpc within ~1 kpc, esti- 
mated from its total stellar mass and assuming that the delensed HST 
1.6 um emission traces the 2. distribution. This is about 1/3 of Z. in the 
densest stellar systems of massive compact elliptical cores inthe local 
Universe“. The cold molecular gas surrounding the centre at 
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Fig. 3 | ERIS spectra of PJ0116-24. a-c, Rest-frame optical lines integrated 

over the entire Einstein ring. a, HB. b, [O 111] AA4959,5007.c, Ha and [N 11] 
AA6548,6584. We fitted the Ha and [N 11] AA6548,6584 lines each with a narrow 
(FWHM = 240 km s”) anda broad component (FWHM = 500 kms”). d, Balmer 
decrement (Ha/Hf line flux ratio) versus stellar mass of PJO116-24 compared with 
SDSS z= 0 galaxies (intensity image with colour indicating the number density), 


MOSDEF z= 1-3 galaxies, and another HyLIRG G165 Arc la and its companion 
galaxy G165 NS 46 at z = 2 (see references and discussion in Methods). 

We show PJ0116-24’s Ha/Hf uncorrected for the stellar Balmer absorption 
(st.uncorr.) as an open pentagon for comparison. Error bars indicate 
observational uncertainties for individual data points (mean values + standard 
deviation (s.d.)). 
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Fig. 4 | Strong-line diagnostics of PJO0116-24’s gas-phase metallicity. Each 
panel indicates a strong-line diagnostic method using the Ha, HB, [N 11], [S11] 
and [O 111] line-integrated flux measurements from the ERIS H- and K-band IFU 
data (Methods). Lines indicate empirical relations for star-forming galaxies from 
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ref. 34 (Methods), with shading indicating the scatter (s.d.) from ref. 34. Solar 
metallicity is indicated by the vertical line. Error bars indicate observational 
uncertainties for individual data points (mean values + s.d.). 


galactocentric radii ~1-2 kpc is abundant and dense, with 
Xeoldgas,1-2kpc % 2 X 10° Mo kpc’, and the SFR surface density there 
iS UseRi—2kpe ~100-200 Mo yr kpc? based onthe delensed distribu- 
tion of dust emission at ~2 mm (with a comparable resolution as CO). 
These 2 cod gas,1-2kpc AN 2 sex 1-2kpc are More than ten times higher than 
those of most local spiral galaxies and are comparable to the maximal 
local gas and SFR surface densities predicted in secularly evolving 
galaxy simulations”. This picture is probably consistent with the forma- 
tion path of compact, quenched bulges. Simulations have predicted a 
‘wet disc compaction’ mechanism“ in which disks fed by cold streams 
and undergoing violent disc instability will drive gas to inspiral towards 
the centre due to dissipative processes. Soon after a high central mass 
surface density is accumulated (X, c1kpe % 10° ° Mo kpc)", internal 
quenching occurs, which consumes the cold gas inside-out and forms 
a massive, compact central bulge**". PJO116-24 appears to have been 


caught in the middle stage of this wet compaction, namely the 
inside-out quenching. This scenario is further supported by the super- 
solar metallicity of this galaxy, given that about half of the stars in the 
Galactic bulge have indeed supersolar metallicity”. 


Asecular evolution pathway for massive HyLIRGs 

PJO116-24 is a unique hyperluminous rotating disk appearing as the 
brightest HyLIRG Einstein ring in the southern sky, and it is one of the 
rare cases studied in detail for several gas phases. Unlike almost all 
other extreme HyLIRGs, which are major mergers’ "4*8, PJO116-24 
does not obviously have massive companions or disturbed kinematics 
as evidence for major mergers. Recently, only one other lensed HyLIRG 
Einstein ring, named 9i09 or PJO20941.3 + 001559 (at z = 2.55, with 
LL p= 1.3 x10 L, and = 10, also in the PASSAGES selection)****° has 
shown evidence of circular rotation, with Umax/Oo,mol. gas % 4-9 + 0.7 
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from the cold-gas tracer CO(4-3) observations at ~360 pc delensed 
scales*’. To date, these are the only two known HyLIRGs with evidence 
of circular rotation over ~10* deg? in the sky. Compared to 9i09 
(ref. 50), PJO116-24’s CO data have a slightly higher spatial resolution 
and reveal much richer substructures, and PJO116-24’s near-IR IFU 
observations probe the heavily attenuated rest-frame optical lines that 
have rarely been imaged in other HyLIRGs. 

The robust confirmation of PJO116-24 as the most rotationally 
supported HyLIRG (U;o¢/00,mol.gas = 9-4 + 0.9 and Vrot/Oo,ion.gas = 
5.7 + 0.6) from this work is key evidence suggesting that secular evolu- 
tion, that is, without recent major mergers, can be responsible for 
maximal star formation in the Universe (intrinsic SFR > 1,000 Mo yr 
in our case PJO116-24 has SFR ~ 1,500 Mo yr~). This scenario has been 
well predicted in simulations before” and is a decisive complement to 
the established paradigm of the evolution of massive (M, > 10'°° Mo, 
gas-rich, turbulent, Toomre-unstable, disk galaxies at z = 1-3 
(refs. 19,26-29) studied at kiloparsec scales. 

In our state-of-the-art three-dimensional (3D) dynamical model- 
ling directly inthe image plane (Methods), we corrected for the lensing, 
spatial beam-smearing, line instrumental broadening and 3D projection 
effects and robustly constrained that PJ0116-24 is genuinely a massive, 
turbulent, rotating disk with baryonic mass Mgisk ayn © ine” Mo 
and intrinsic cold-molecular and warm-ionized turbulence dispersion 
Oo,mol.gas ¥ 43 + Skms~'and 09 jon. gas % 69 + 8 km s™'. Moreover, it is 
one of the rarest HyLIRG Einstein rings in the Universe. This study 
presents a unique case that not only exemplifies the secular evolution 
pathway for HyLIRGs but also serves as the most detailed 
high-resolution, multi-tracer study of the HyLIRG population, whichis 
nowadays routinely found in existing and future wide-area surveys. 


Methods 

VLT ERIS IFU Science Verification Program 

Through the ERIS Science Verification Program (ID 110.258 S, PID. Liu), 
we observed PJ0116-24 in the H- and K-bands in seeing-limited condi- 
tions with the high-spectral-resolution H_short and K_short gratings 
and the 250 mas pixel scale. The two gratings have a spectral resolution 
of R= 11,200 and cover A = 1.46-1.67 um (1.93-2.22 um) in the H-band 
(K-band). The observations were taken on 2 December 2022 during the 
end of ERIS commissioning (https://www.eso.org/sci/activities/vltsv/ 
erissv.html), as for the other ERIS Science Verification Programs, which 
aimed to verify the capabilities of this new instrument. 

The observation blocks included a first acquisition to a 
nearby bright ‘cal_PSF’ star (Gaia DR3 5040854651181560064) 
and then a sequence of object exposures (J2000 right ascension 
01h 16 min 46.7792 s and declination -24° 37' 02.518”) and sky expo- 
sures (J2000 right ascension 01h 16 min 49.1979 s and declination 
-24° 37’ 24.918”) with small dithering offsets (+0.5”) to allow for a 
better combined spatial sampling and to avoid the effects of bad or 
hot detector pixels. For each H- and K-band observation, we obtained 
eight object exposures and four sky exposures with a sequence of 
object-sky-object and object-sky-object dithering. The dither pat- 
tern of the object exposures was [(—0.5”, -0.5”), (+0.0”, —0.5”), (+0.5”, 
-0.5”), (+0.5”, +0.0”), (+0.5”, +0.5”), (+0.0”, +0.5”), (-0.5”, +0.5”), (-0.5”, 
+0.0”)]. The sky exposures pointed towards the same blank sky. Each 
exposure had 300 s of Detector Integration Time (DIT). 

The latest ERIS pipeline was used to perform flat fielding, distor- 
tion correction, wavelength axis calibration, atmospheric OH line 
fitting and subtraction, sky exposure subtraction from object expo- 
sure data, and correction for differential atmospheric refraction. We 
calibrated the astronomical data unit to physical flux scale conversion 
based on ‘telluric’ stars observed with the same filter, spectral resolu- 
tion and pixel scale as our observations. We collected all available 
telluric star observations from April 2022 to July 2023 and determined 
a mean physical flux scale conversion spectrum by normalizing to a 
telluric star’s theoretical black body (with temperature from stellar 


type) or template SED (see the European Southern Observatory (ESO) 
website; for example, https://www.eso.org/sci/facilities/paranal/ 
decommissioned/isaac/tools/lib.html). The variation of the flux scale 
conversion spectrum was ~20%-50%, which could be due to the vary- 
ing atmospheric conditions of observations not considered for inthe 
pipeline (the conditions considered were air mass and DIT but not the 
seeing or water vapour). 

After obtaining individual calibrated 300 s exposure data cubes, 
we combined each cube witha custom script that performs sigma clip- 
ping for bad pixels, image reprojection and weighted-averaging chan- 
nel by channel. The weighting was based on the spectral pixel (spaxel) 
errors determined from sigma-clipped statistics. The final combined 
image cube was smoothed by a 2 pixel boxcar kernel to increase the 
S/N. The PSF was determined by fitting a 2D Gaussian to the ‘cal_PSF’ 
star, which was also smoothed by a 2 pixel kernel. The PSF FWHM was 
0.7” in the K-band data and 1.2” in the H-band data. 

Emission line spectra were extracted within a polygon region that 
corresponds to the full Einstein-ring area seen in the HST stellar emis- 
sion when smeared to the low resolution (Supplementary Fig. 1). The 
uncertainties for the spectra were extracted correspondingly using the 
aforementioned error cube, as shown in red in Supplementary Fig. 1 
(right). Spikes indicate channels contaminated by the atmospheric OH 
line (for example, ref. 51). Then the one-dimensional flux-calibrated 
spectra were fitted simultaneously with lines and a polynomial con- 
tinuum, for each H- and K-band. The line flux errors were determined 
based on the fitting method, for which we used the astropy.modeling 
software package” ™ with a fitting type of LevMarLSQFitter. 

The Ha and [N 11] AA6548,6584 lines were also fitted with a broad- 
line and a narrow-line component for each of the lines, as shown in 
Fig. 3b, due to their high S/N. We assumed that the broad-line com- 
ponents have the same linewidth and line-centre relative offsets, 
whereas the narrow-line components have a different linewidth and 
set of line-centre offsets. We found redshifted, broad-line components 
for all Ha and [N 11] AA6548,6584 data points. The broad-line FWHM 
was -510 km s”, about twice of that of the narrow-line FWHM, which 
was ~240 km s™. Instead of considering it as evidence of an outflow, we 
verified that it was mostly due to the asymmetric lensing magnification 
of the blueshifted and redshifted parts of the galaxy, as described in 
the main text. In our later rotating-disk forward modelling, we were 
able to recover avery similar asymmetric line profile in the image plane 
froma modelled pure Sérsic-profile rotating disk in the source plane. 
Therefore, the global asymmetric spectrum does not provide con- 
vincing information of an inflow or outflow. A future high-resolution 
near-IRIFU follow-up will shed light on any inflow or outflow signal for 
the ionized gas across the galaxy and at the AGN. 

Then, we corrected the stellar Balmer absorption for the Ha and 
HB line fluxes using our SED models based on BCO3, Bruzual and Char- 
lot’s stellar population synthesis code (‘Fitting the SED’)**. We used 
the high-resolution stellar library to generate a model SED with stellar 
Balmer absorption and measured the Ha and Hf stellar absorption 
equivalent widths as 4.1and 5.4 A, respectively. These equivalent widths 
were multiplied with an emission filling fraction of 0.36 (0.23) for Ha 
(HB)*, multiplied with the corresponding stellar continuum flux densi- 
ties and then added to the directly measured line fluxes. 


ALMA observations 

We obtained several ALMA observing programmes (ID 2017.1.01214.S, 
PIM. Yun; ID 2019.1.01197.S, PI P. Kamieneski; and ID 2021.1.00353.S, 
PI K. Harrington) to observe the molecular gas CO line emission and 
dust continuum in the PASSAGES strongly lensed HyLIRG sample. 
The CO(3-2) data analysed in this study was combined with data from 
programme 2017.1.01214.S with 7.08 min on-source integration at 
0.31” angular resolution and with data from programme 2019.1.01197.S 
with 4.56 min on-source integration at 0.35” resolution and 19.08 min 
on-source integration at 0.10” resolution. 
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In this work, we calibrated the raw ALMA measurement sets with 
the standard calibration pipeline with the Common Astronomy Soft- 
ware Applications (CASA). We then regridded the visibility data toa 
common spectral axis using the CASA mstransform task with a chan- 
nel width of 30 kms ‘and subtracted the continuum with the CASA 
uvcontsub task. The CO(3-2) image cube was produced by running 
the CASA tclean task, first with multiscale models with a threshold 
downto 3gand then by adding hogbom (single-scale) models downto 
lo, where gis the root mean square of the cube determined after a few 
major cycle iterations. A clean mask was applied based onthe smoothed 
version of a first-run CO(3-2) image cube. This tclean approach has 
been found very useful for imaging galaxies with extended emission 
structures*®. We used the natural weighting, which gives the best S/Ns 
at a sufficiently small, synthesized beam (FWHM = 0.19” x 0.16”) for 
our analysis. The maximum recoverable scale was 7.5”, which was suf- 
ficient for recovering the extended emission from PJO116-24. We also 
produced Briggs-weighting clean cubes, which have higher angular 
resolutions and better suppressed sidelobes but are noisier. Thus, we 
used them mainly for visually inspecting the emission peaks in our 
lensing modelling. 

As our CO dataset includes several array configurations, the syn- 
thesized clean beam deviates from a 2D Gaussian shape. There is a 
factor of 2-3 difference between the effective area of the clean beam 
and the dirty beam, leading to an incorrect flux scale in the clean resid- 
ual (that is, emission that was not contributed from any clean model 
component) and affecting the total flux measurements by a few tens 
of percent depending on the clean depth and the clean-to-dirty beam 
area ratio £. Therefore, a correction was needed to scale down the flux 
within the clean residual by the factor £, whichis called the JvM correc- 
tion”. We applied this correction to all our imaged cubes by computing 
Ieteanimage — (1 — £) Ictean residua Where € is measured from the ratio between 
the clean beam area and the area of the dirty beam image integrated 
from the centre to the first null radius. 

The ALMA dust continuum images at various bands were produced 
similarly and used to constrain the SED of PJO116-24 (‘Fitting the SED’). 
We compared the dust continuum images and fluxes made by us and 
made by previous studies*ć and found good agreement. 


VLT MUSE observations confirmed the lens redshift 

Next we report VLT MUSE IFU observations around PJ0116-24, which 
spectroscopically confirmed the foreground lens to be a massive 
early-type galaxy (J2000 right ascension 01h 16 min 46.7963 s and 
declination -24° 37’ 02.234”) at z= 0.554. The MUSE observations were 
obtained in a filler programme (programme code 110.23, PI F. Bian) 
under poor seeing weather conditions during 5 to 6 October 2022, 
with a seeing of 0.86”-1.80”. The MUSE data cover a wavelength range 
of 4,800-9,300 A with a spectral resolution of R = 1,770-3,590. The 
wide-field mode offers a field of view of 1’ x 1’ witha pixel scale of 0.2”. 
The data reduction used the standard MUSE data reduction pipeline 
provided by ESO (https://www.eso.org/sci/software/pipelines/; see B. 
Alcalde Pampliega et al., manuscript in preparation). 

Supplementary Fig. 2 shows the MUSE spectrum extracted at this 
foreground lens within an aperture of 1.3” and a fitted spectrum using 
the pyplatefit Python package’s fit_spec function® and with identified 
line names labelled. Strong absorption lines from old stellar popula- 
tions are clearly visible in the spectrum, especially the prominent Ca 
doublet. Thus, we determined the redshift of the massive lens galaxy 
to be z= 0.554. A few more sources have a MUSE spectrum indicating 
that they are at z= 0.554. Thus, they could be inthe same galaxy group 
or cluster as the massive galaxy acting as the primary lens of PJO116-24 
(B. Alcalde Pampliega et al., manuscript in preparation). However, the 
closest second member of this potential galaxy group or cluster is at 
8.2” away from the massive lens galaxy, so that they have little impact 
onthe lens modelling for PJO116-24. We considered a dark matter halo 
plus an external shear in addition to the massive lens galaxy itself (see 


the next section). Thus, we took the effect of this potential galaxy group 
or cluster into account. 


HST and Gemini optical and near-IR data 

The HST WFC3 F160W data are from the observation programme 
GO-14653 led by PIJ. Lowenthal (J. Lowenthal et al., manuscript in 
preparation) and ref. 6 has details of these observations. The image 
was dithered with five positions and combined to achieve a pixel size of 
0.065”. The PSF is ~0.15”. The point-source So depth is my, = 27.4. There 
are no other HST or James Webb Space Telescope (JWST) observations 
for PJO116-24 as at the time of writing. 

For PJ0116-24, there is also Gemini multi-object spectro- 
graph r- and z’-band imaging from the Gemini south programme 
GS-2018A-Q-216 (PIJ. Lowenthal). The Gemini imaging for the whole 
PASSAGES sample will be presented in O. Cooper et al. (manuscript in 
preparation). A summary has also been presented in ref. 6. The images 
have a pixel size of 0.16” and the seeing-limited PSF is ~0.7”-0.8”. 


Lens modelling 

PJO116-24 is primarily lensed by a foreground elliptical galaxy and its 
dark matter halo. We used the following steps to construct our lens 
model. 


(1) We corrected the HST astrometry using the Gaia DR3 star 
catalogue as a reference and then performed source fitting of 
the HST images. We manually added prior sources at HST emis- 
sion peaks and used GALFIT for the prior source fitting with 
Sérsic models. The prior sources include both the foreground 
elliptical galaxy and the lensed PJO116-24. As the foreground 
elliptical galaxy has strong and extended HST emission, we 
used two Sérsic models to model it best. Once we had obtained 
the best-fitting parameters, we ran GALFIT again but only to 
model the foreground elliptical galaxy without fitting to obtain 
the foreground-subtracted HST image only with PJO116-24’s 
emission. The foreground-subtracted image is shown in Fig. 1 
and Supplementary Fig. 3 (see also refs. 5,6, and J. Lowenthal 
et al., manuscript in preparation where the original HST image 
is shown). 

(2) We identified emission peaks (‘knots’) in HST and ALMA channel 
CO(3-2) maps and used their sky coordinates as the input for 
our later lens model optimization. As demonstrated in ref. 61, 
which analysed a x10 less luminous, massive, strongly lensed 
galaxy JO901-1814 in a similar way, it is very important to have 
both HST and ALMA knots as they provide highly complemen- 
tary spatial coverage. The ALMA channel maps are critical be- 
cause they provide unique velocity information for the pairing 
of knots. Supplementary Fig. 3 shows the knots identified for 
HST and CO emission peaks as well as the CO and dust intensity 
and CO velocity maps. 

(3) With the knot positions as the input, we used the glafic soft- 
ware? to optimize the lens model parameters. Our model 
set included a dark matter halo with a Navarro—Frenk-White 
profile, the foreground elliptical galaxy as a singular isother- 
mal ellipsoid model and external shear to allow for better op- 
timization (see also ref. 61). We ran the glafic default minimum 
chi-square algorithm to determine the best-fitting lens model 
parameters. The glafic input files and lens model parameters are 
publicly available with this work. 

(4) After obtaining the best-fitting lens model parameters with the 
minimum chi-square algorithm, we used a Markov chain Monte 
Carlo (MCMC) algorithm to explore the parameter space and 
to verify that the parameters were tightly constrained, follow- 
ing ref. 61. The MCMC fitting gives the probability distribution 
of each lens model parameter and verifies that the minimum 
chi-square lens model is within the 1o confidence regions. The 
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MCMC-based fitting results are also publicly available with 

this work. As discussed in ref. 61, MCMC fitting is useful for 
obtaining the probabilities and uncertainties, but the minimum 
chi-square lens model provides the best fits to the knots. There- 
fore, we used the minimum chi-square lens model as our best 
model for the further analyses. We visually verified the recon- 
structed source-plane maps to ensure the reliability of our best 
model (as can be seen in Figs. 1 and 2 for the different tracers). 
In addition, we also compared our delensed maps with those 
from Kamieneski et al.° using the lenstool software and based 
onthe HST and ALMA 0.4” resolution 1.1mm dust continuum 
maps. We found very similar delensed maps. Kamieneski et al.° 
reported a magnification of 44 imm = 7 for the 0.4” resolution dust 
continuum map, and their model also implies Hysr = 13.8 for the 
same 0.15” resolution HST map as used in this work (private 
communication). Their Hust is, thus, within 13% of the Hust = 15.9 
in this work, which was based on the 0.15”-0.19” resolution HST 
and ALMA CO channel maps. 


After obtaining the best lens model, we converted the image-plane 
data into the source plane, namely delensing. We used the mesh grid 
from the best-fitting lens model, which describes how rectangle foot- 
prints (‘cells’) in the image plane are mapped into polygons in the 
source plane. We looped over mesh grid cells and co-added the emis- 
sion from the image-plane pixels (in each cell’s rectangle footprint) 
into pixels in the source plane (in each cell’s polygon footprint) with 
linear interpolation inside each cell. For the CO and Ha cubes, we 
always first reconstructed the source-plane data cubes then extracted 
the line amplitude, velocity and velocity dispersion properties with 
one-dimensional Gaussian profile fitting for each spaxel in the cube, 
with either the scipy.optimize.curve_fit or pymc3 Python package. 
Note that the source-plane images were mostly for visualization and 
were not directly used for our kinematic modelling. Instead, the lens 
model’s mesh grid was used by our kinematic modelling software 
DysmalPy+Lensing to allow direct fitting of the image-plane data 
(‘Dynamical modelling with DysmalPy+Lensing’). 

As previously mentioned, we measured the ‘delensed beam’, thatis, 
the angular resolution in the source plane, by generating 2D Gaussian 
beams in the image plane at various positions then delensing each of 
them into the source plane and analysing their shapes. For image-plane 
2D Gaussian beams (FWHM = 0.16”) at positions Aland A2, we measured 
adelensed beam FWHM of 0.088” x 0.055” and 0.149” x 0.039”, respec- 
tively. Their minor-axis FWHM, thus, correspond toa physical scale of 
450 and 320 pc, respectively, atz = 2.125 (8.18 kpc arcsec”). Close to the 
critical curves, the best delensed angular resolution can reach <100 pc 
along the northern CO-bright arc, for instance, along the C2-C1-B4 
knots in Supplementary Fig. 3. This happens to be the location of the 
CO arm-like substructures, which allowed us to spatially resolve these 
disk substructures even better than the centre, as shown in Fig. 1d. 

In Fig. 1d, we also see some other Ha blobs having no HST or CO 
counterparts. The brightest oneis north of clump B, as mentioned in the 
caption. We made tests to verify that this Ha blob is probably an artefact 
of delensing. The relatively poor angular resolution of our ERIS data 
created beam-smeared emission and noise fluctuations. Those that are 
scattered right across the critical curves, for example, around B2, were 
mapped to the exact position of the artificial Ha blob. Higher angular 
resolution ERIS observations with adaptive optics will probably solve this 
issue, as beam-smearing across critical curves will be largely reduced. 


Comparison with Balmer decrements in the literature 
Figure 3 compare the Balmer decrement of PJO116-24 with the follow- 
ing galaxy samples: 


e local SDSS galaxies from the SDSS DR8 MPA-JHU catalogue”, 
with type ‘GALAXY’ and with S/N = 10 in both Ha and Hf line fluxes 


e MOSDEF high-z galaxies from the MOSDEF survey®*, with S/N > 4 
in both Ha and HB line fluxes 

e G165 Arc la, a HyLIRG that is also in the PASSAGES sample, and 
its companion galaxy NS 46; they are both in the G165 field” 


There are also Balmer decrement measurements for several hun- 
dred galaxies of different types, for example, local ULIRGs” and 
star-forming galaxies atz = 0.75-7 (refs. 71-76). Local ULIRGs are rare, 
merger-driven starbursts. A few known ULIRGs (Arp220, UGC05101 
and IRAS12112 + 0305) have an extreme Ha/Hf = 9-20 (ref. 70). Never- 
theless, the majority of high-z galaxies with rest-frame optical line 
measurements do not havea global Ha/Hf as high as that of PJO116-24, 
and their stellar masses and SFRs are, in general, much lower. An excep- 
tional case was reported recently by Frye et al.°’, G165 Arc 1a, a HyLIRG. 
This HyLIRG is also in the PASSAGES sample and has JWST NIRSpec 
multi-object spectroscopy indicating that Ha/Hf = 11.7 + 3.1 (corrected 
for stellar Balmer absorption). It is probably less massive (intrinsic 
stellar mass M, ~ 10'°° Mg) and much less magnified (4 ~ 2.7) than 
PJO116-24. It also appears to be interacting with a companion galaxy 
named NS 46, which is 1,400 km s” offset from the systematic velocity 
of G165 Arc la and is an even less-massive galaxy (intrinsic stellar mass 
of -10° M,,; ref. 69). Both G165 Arc 1a and PJ0116-24 are HyLIRGs with 
extremely high Ha/Hf, but their masses and merger history are very 
different. Thus, both are unique cases for understanding extreme dust 
attenuation. 


Ionization source and gas-phase metallicity diagnostics 

With the detected nebular emission lines in PJO116-24, we analysed 
the ionization source using a BPT diagram”. Supplementary Fig. 4 
shows the [N 11]-BPT” and [S 11]-BPT” diagnostics. For the [N 11]-BPT, 
we adopted the starburst-AGN classification curve from ref. 79, 
the star-forming and composite star-forming division curve from 
ref. 80 and the separation curve between Seyfert and LINER-type AGNs 
from ref. 81 (see also ref. 82). For the [S 11]-BPT, we used the ref. 82 
classifications. 

Our galaxy falls into the star-forming/composite regime. That is, 
the emission lines are mostly powered by star-forming H II regions at 
global scales. This diagnostic is conservative, as refs. 83,84 and many 
studies have found that the classification lines for high-redshift galaxies 
need to move towards higher [N 11]/Ha and [O 111]/HB, which puts our 
galaxy further inside the star-forming regime. 


Interstellar medium mass from CO and [C1] line modelling 

The mass of the interstellar medium (ISM), whichis considered to be 
molecular gas in this study, was measured for PJO116-24 as 
log(uMsy/Mo) = 12.4970 65 from an updated MCMC run of the radiative 
transfer modelling in Harrington et al.* (to be presented in K. Har- 
rington et al., manuscript in preparation). We re-imaged all available 
ALMA CO and [C 1] line data using the cleaning and JvM correction 
method as mentioned above in ‘ALMA observations’. This included the 
highest-resolution CO(3-2) data used in this study and lower-resolution 
CO(4-3), CO(7-6), CO(9-8), [C 1](1-O) and [C 1](2-1) data observed by 
ALMA programmes 2017.1.01214.S (PI M. Yun), 2019.1.01197.S (PI P. 
Kamieneski), 2021.1.00353.S (PI K. Harrington) and 2021.2.00062.S (PI 
D. Riechers). The single-dish CO(1-0), CO(3-2), CO(6-5), CO(7-6), 
CO(8-7), CO(9-8) and [C 1](2-1) line fluxes reported in Harrington 
et al.” were also used, and these were consistent with our interfero- 
metric measurements when both exist. 

The Harrington et al.” radiative transfer modelling assumes a 
turbulent gas density probability distribution function to describe the 
ISM. CO and [C1] molecular line emission was modelled by integrating 
the contributions from different gas densities. The gas kinetic tem- 
perature (Tn) was modelled as a function of the H, density (m,,). The 
parameters of the temperature-density function, CO and [C 1] abun- 
dances, ny, cloud radius properties and turbulence velocity dispersion 
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AV up Were all free parameters. Harrington et al.* reported 


log(uMism/Mo) = 12.67 + 1.15 for PJO116-24, witha relatively large uncer- 
tainty mostly from the limited CO and [C 1] lines and the degeneracy 
among parameters in that study. Our updated modelling included more 
data so that the resulting M,.y has a much tighter constraint. 

There are still a few types of uncertainty that our modelling did 
not account for. For instance, the assumed monotonic form of the 
temperature-density function is probably too simple to account for 
the real ISM conditions, as there could bea range of gas temperatures 
for a given gas density. Indeed, by running a simpler two-component 
radiative transfer model following Liu et al.%, we found a slightly 
smaller Mısm (log(uMism/Mo) © 12.38 + 0.15 whichis at the lower lobound- 
ary of the updated Harrington et al. radiative transfer modelling) due 
toawarmand very dense gas component in addition to a cold and dense 
component. However, the two-component modelling achieves poorer 
fitting to the observed data, especially at mid-/CO lines. Thus, we used 
the updated Harrington et al.® radiative transfer modelling result. 


The scenario of a giant star-forming clump 

Clump B, whichis labelled in Fig. 1, was detected in the HST and Gemini 
images, as well as inthe ERIS Ha map and CO channel map but notin the 
dust continuum images. It is quadruply lensed into the B1-B4 counter- 
images. As discussed in Kamieneski et al.°, this lensing configuration 
indicates that the clump should be at the same redshift as the Einstein 
ring. The Ha line moment-O map (Fig. 1d) clearly confirms this clump 
as all four counterimages were detected, and it indicates that it is form- 
ing stars. The CO(3-2) emission was also detected at a significance of 
~2-4o in the spectra extracted from the clump at the B1-B4 locations, 
as shown in Supplementary Fig. 6. The integrated CO(3-2) fluxes of the 
clump at the B1-B4 locations were 0.19 + 0.10, 0.29 + 0.10, 0.40 + 0.16 
and 0.23 + 0.06 Jy kms”, respectively. They are in line with the mag- 
nification factors estimated for the B1-B4 locations (based on our 
best-fitting lens model magnification map smoothed to the CO(3-2) 
angular resolution), which are ~9-10 at B1 and B4 and ~14-15 at B2 
and B3. Given that the total integrated CO(3-2) flux of PJO116-24 was 
66.4 + 8.5 Jy km s” from our new data reduction (which is consistent 
within ~1o to the measurement 51.10 + 10.22 Jy km s” reported in Har- 
ringtonetal.*), the clump’s total CO emission is only -1.6% of the total 
galaxy’s CO emission. If the CO emitting gas in the clump and across the 
PJO116-24 disk has similar excitation conditions and, thus, a CO-to-H, 
conversion factor, then the CO line flux ratio is also a molecular gas 
mass ratio, indicating that the cold ISM mass is negligible in the clump 
compared to the entire galaxy. 

The HST 1.6 um emission of the clump is about 13% (1/7) of the 
whole PJ0116-24, but this is a rest-frame 5,000 A light ratio instead of a 
mass ratio. Compared to the central kiloparsec of PJO116-24, the clump 
is ~1 mag bluer inthe z- H colour (roughly based on the Gemini z-band 
and HST H-band images, although their angular resolutions are very dif- 
ferent). Based onthe BCO3 (ref. 55) stellar population synthesis models, 
~I mag bluer z - H colour corresponds to a H-band mass-to-light ratio 
lower by a factor of 5 under solar metallicity (for example, based onthe 
SED templates in Liu et al.°’; see also the appendix in ref. 88). Thus, the 
clump probably has a stellar mass ratio of less than 1-3%, which also is 
negligible compared to the entire galaxy. 

The [N 11] line map does not show detections at the four clump 
counterimages as clear as in the Ha line map, probably due to the 
relative line S/N. The [N I1]/Ha ratio could be useful for probing the 
metallicity of the clump, but realizing an accurate ratio is subject toa 
future high-resolution IFU follow-up study. 

Given the aforementioned environment of PJ0116-24 and the 
detection of stellar light and cold and ionized gas in the clump, we 
believe that the origin of the clump is very probably in situ. It may 
have formed inside the massive turbulent gas disk due to gravita- 
tional instability and was fed by inflowing gas from the dark matter 
halo?30318°102 Such Ha-bright, giant, star-forming clumps have been 


found in many massive unlensed galaxies at z = 1-3, as studied with 
adaptive-optics-assisted near-IR IFU observations”? ?°%04 (see also 
the review by Förster Schreiber and Wuyts’’). An ex situ scenario is 
less probable, as it requires the clump to be a low-mass (~10° M,) blue 
starburst with a high SFR but little dust while having a relative motion 
that is co-rotating with the disk of PJ0116-24. 


Dynamical modelling with DysmalPy+Lensing 

We performed 3D forward kinematic modelling that incorporates 
the lensing effect using the Dysmal software (refs. 19,26-28, 
105-107,108-111) and the recent DysmalPy version’” (see also 
refs. 61,104,113-118) with the Lensing module developed by ref. 61. The 
advantage of this forward-modelling approachis that the beam-smearing 
and lensing are simultaneously taken into account. This avoids the 
complication of modelling a system to fit the delensed data that have 
a spatially varying PSF or inaccuracies resulting from approximating 
the variable resolution with a uniform value. We briefly describe the 
DysmalPy+Lensing forward modelling below (see ref. 61 for more details). 

The source-plane galaxy model is a combination of two geomet- 
rically thick Sérsic-profile components, namely a disk and a bulge, 
and a dark matter halo with a Navarro-Frenk-White profile. The disk 
and bulge components represent the baryonic matter in the galaxy 
and have a systematic rotation and an intrinsic velocity dispersion. 
These two Sérsic-profile components may not physically represent 
a clearly distinct disk and bulge as in local galaxies, but this allows 
for a wide flexibility of the mass distribution model and, hence, the 
shape of the rotation curve. The 3D geometry of the disk + bulge is 
defined by an inclination angle (i) and a position angle (PA), which 
we determined from the delensed HST image and CO velocity field. 
The disk and bulge half-mass radii Re aisk and Re vuige Were initialized 
with values obtained from fitting the delensed HST + CO mass radial 
distribution with two Sérsic profiles (see, for example, ref. 61). Here, 
Re disk ~ 3 kpc and Re vuige ~ 0.6 kpc. The disk size is at the lower boundary 
of the size-mass relation of main-sequence star-forming galaxies” 
and is in the regime for the transition to compact quiescent galaxies 
of the same mass. The disk has a scale height that is fixed to 1/4 of its 
scale length, whichis atypical value for a thick disk atz = 2 (refs. 112,117). 
The bulge-profile mass component was found to contribute little to 
the kinematically fitted baryonic mass as the fitted bulge-to-disk ratio 
(B/T) was less than 0.1. 

We performed the DysmalPy+Lensing fitting for a spatially and 
spectrally smoothed (0.4”) ALMA CO(3-2) data cube to increase the 
sensitivity to the extended, outer rotation curve. Supplementary Fig. 7 
shows the 2D projected velocity maps of the data and the best-fitting 
model, which agree well within the uncertainty. The velocity disper- 
sion in the data was generally fitted well in the disk but less well near 
the galaxy centre, where our idealized model does not account for 
any noncircular motion, which could increase the dispersion and also 
slightly twist the velocity field’. 

Our fitting method used MCMC to sample the high-dimensional 
parameter space. This approach allowed the disk + bulge mass Mask, ayn» 
Re, B/T, intrinsic velocity dispersion do, the dark matter halo virial 
mass Mpm, ayn, Land PA to vary (as well as allowing the geometric centre 
in the source plane to shift slightly). We applied a Gaussian prior on 
log(Maisk, dyn) (with a mean of -11.5 and s.d. of -0.35) guided by the pho- 
tometrically measured ISM and stellar masses. We found that using a 
flat prior resulted in a slightly lower log(Mgisx, ayn) Within 10 error. We 
also applied Gaussian priors on siniand PA based on photometrically 
determined values. 

The DysmalPy+Lensing fitting algorithm generated a 3D galaxy 
model and circular velocity per voxel in the source plane, which was 
projected onto the sky and deflected to the image plane. Then, it gen- 
erated the effects of the PSF and the line spread function. It extracted 
the 2D velocity and dispersion fields, as for real data, and finally, it 
calculated the differences between the synthetic model data and the 
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real data (velocity and dispersion fields). We ran the MCMC fitting with 
250 walkers, each with 300 steps after about 30 burn-in steps. Supple- 
mentary Fig. 8 shows the MCMC posterior probability distributions of 
our fitting, with a few key parameters highlighted. All parameters were 
reasonably well constrained, and the residuals could not be further 
improved given our idealized mass model assumptions (introducing 
noncircular motion could bea future direction to explain the residuals 
seen in Supplementary Fig. 7; see ref. 104). 

Our obtained disk baryonic massis Maisk,dyn © 10 Mo whichis 
not far from thesum of the ISM andstellar masses (Mism + M.) ~ 10" Mo 
where Misy ~ 10"7°*°? Mand M, ~ 10""°*°? M, Theexcess of (Misy + M.) 
compared to Maisk, ayn İS at $ 10 significance. This could be due to 
several assumptions or uncertainties unaccounted for, such as the 
star-formation history in our spectral energy fitting, the temperature- 
density coupling function in our ISM mass modelling, the initial mass 
function (IMF), which may change in high-z starbursts”°, and the con- 
tribution ofthe AGN tothe continuum and line emission. All these uncer- 
tainties may add up to -0.5 dex for the masses, so that the marginal 
difference between Mäisk, ayn and (Mism + M.) is not a significant issue in 
this study. 


11.2940.3 


Fitting the SED 

We ran multi-wavelength SED fitting including optical, IR, (sub) mil- 
limetre and radio bands to constrain the stellar and dust properties of 
the galaxy using the MiChi* code”. The photometry data were mostly 
from refs. 4-6, except that we reduced the ALMA 180, 260 and 395 GHz 
dust continuum consistently (Supplementary Fig. 3). 

The SED model is a composition of a stellar component generated 
using the BCO3 code” witha constant star-formation history and solar 
metallicity, a mid-IR AGN using empirical templates from ref. 121, a 
Draine and Li dust model heated by a physically driven interstellar radia- 
tion field’ (updated in ref. 123) and a simple power-law synchrotron 
radio component. We assumed a Chabrier IMF‘. A bottom-heavier 
Salpeter IMF'” would result in a stellar mass and SFR a factor of 1.7 
larger. A top-heavier IMF would result in a smaller stellar mass and 
SFR. The MiChi’ code” fitted the four components simultaneously and 
explored the multi-dimensional parameter space with Monte Carlo 
realizations to sample the chi-square-based probability distribution 
function of each parameter. 

Supplementary Fig. 10 shows the SED of PJO116-24 and the probabil- 
ity distributions of the apparent (magnified) stellar mass, stellar con- 
tinuum dust attenuation E(B - V).,,., IR luminosity and dust mass 
(Table 1). Adopting a usual stellar continuum to nebular dust attenuation 
ratio Nne» = 0.5 (refs. 126-129), the SED-based E(B — V) ney, sep = 1.00 + 0.26 
which is fully consistent with the Balmer-decrement-based measure- 
ment reported inthe main text. 


Data availability 

ESO VLT ERIS and MUSE data are publicly available at the ESO archive 
(https://archive.eso.org/cms/eso-data/instrument-specific-query- 
forms.html). ALMA data are publicly available at the ALMA Science 
Archive (https://almascience.eso.org/aq/ with the identifiers ADS/JAO. 
ALMA#2017.1.01214.S and ADS/JAO.ALMA#2019.1.01197.S). HST data 
are available at the Mikulski Archive for Space Telescopes (https://mast. 
stsci.edu/). Our ALMA CO and ERIS Ha + [N 11] cubes and the glafic lens- 
ing modelling files are available via figshare at https://doi.org/10.6084/ 
m9.figshare.25359613 (ref. 130). 


Code availability 

ERIS and MUSE data reduction pipelines are available at the ESO VLT 
Instrument Pipelines webpage (https://www.eso.org/sci/software/ 
pipelines/). Custom scripts to combine the single-exposure ERIS data 
cubes are available upon request. The ALMA data reduction software 
CASA” is available at the CASA download webpage (https://casa.nrao. 
edu/casa_obtaining.shtml). Custom scripts to combine the visibilities 


and perform imaging are available upon request. Our lens modelling 
glafic input and best-fitting files are made publicly available with this 
paper. The Dysmal kinematic fitting tool (with its Lensing module) 
will soon be made public™ and is also available upon request before 
its public access. Other public tools used for the data analysis are 
astropy” “, astroquery'”, GALFIT’**, glafic®’, matplotlib’’°, MiChi? 
(refs. 87,137), numpy’’ and scipy’”. 
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